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1. Introduction 

The pions are supposed to be the Nambu-Goldstone bosons associated with the spontaneous 
chiral symmetry breaking in QCD. The three pions {n + , 7T°, tt~) form an isospin triplet of flavor 
SU(2) symmetry. Among the three pions, 71° is most unstable one, with a lifetime ~ 10~ 9 times 
shorter than that of the other two. Experiments show that the neutral pion decay is mainly an 
electromagnetic process, with the bulk of the decay rate going to two photons. 

At the leading order of QED the 71° — > yy decay width can be expressed as 

7LGC in 

r ^n = -^^(mlAO) , (1.1) 

where a e is the fine structure constant, m % is the neutral pion mass and ^ n 0yy{m 2 7t ,p\,p\) is the 
71° — > yy transition form factor with p\2 the photon momenta. In 1967 Sutherland and Veltman 
argued that in the chiral limit ,^0^(0,0,0) vanishes and the pion can not decay into photons [1,2]. 
At the physical pion mass the predicted decay width is ~1000 times smaller than the experimental 
measurements. Later it was realized that the PCAC relation used in Sutherland and Veltman's 
analysis is only valid at the classical level. In quantum field theory, the conservation of the axial 
current is violated by quantum fluctuations. At the one-loop level, a fermionic triangle diagram 
contributes an extra anomaly term (ABJ anomaly) to the PCAC relation and makes ^^.o 7r (m|,0,0) 
non-zero in the chiral limit [3,4] 

&f y \ = ^(0,0,0) = ■ (1-2) 

In Eq. (1.2) N c is the number of QCD colors and Fq is the pion decay constant F n in the chiral limit. 
It is proved in perturbation theory that higher-loop diagrams do not contribute to J^o^ (0,0,0) [5]. 
As a consequence, the ABJ anomaly gives a rather precise predication for the n° — > yy decay rate. 

In the chiral and on-shell photon limit, the pion decay amplitude is constrained by the ABJ 
anomaly. Away from these limits some corrections from QCD are expected. The recent PrimEx 
experiment at JLab has measured the neutral pion decay width to an accuracy of 2.8% [6]. The next 
stage of this experiment is to achieve a precision of 1.4%. At this level of accuracy the correction 
due to finite quark mass becomes relevant. In this paper we report a model-independent calculation 
of the 71° — > yy form factor and decay width using lattice QCD. The first motivation of this work 
is to determine the finite quark mass correction from first principles. Our second motivation comes 
from the hadronic light-by-light (HLbL) scattering, which is responsible for the second largest 
theoretical error in the determination of the muon g-2. While the direct QCD calculation of HLbL 
is very demanding as it involves a four- vector-current correlation function, the 7T° — > yy form factor 
can be used to estimate the dominant pion-exchange contribution to the HLbL. Thus our calculation 
serves as an intermediate step towards the precise determination of HLbL. 

2. Chiral anomaly on the lattice 

The chiral anomaly is of central importance for the neutral pion decay. When we perform 
a lattice calculation, a natural question is how the chiral anomaly is achieved on the lattice. The 
answer depends on the formulation of the fermion action. In the case of Wilson fermion, the chiral 
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symmetry is explicitly violated and the anomaly is recovered by taking the continuum limit of 
the chiral symmetry breaking term [7]. For the Ginsparg -Wilson fermions, the chiral symmetry is 
preserved in a modified form [8]. The chiral anomaly is introduced by the Jacobian of the chiral 
transformation in this case. When the background gauge field is sufficiently smooth, the Jacobian 
yields the correct chiral anomaly up to discretization effects [9]. 

In our calculation we use the overlap fermion formulation, which is a realization of the Ginsparg- 
Wilson fermion on the lattice. At practically used lattice spacings (~0.1 fm) the gauge field is far 
from smooth and the chiral anomaly may not be guaranteed. Therefore it is important to check 
whether the chiral anomaly is correctly reproduced in our calculation. 

3. Treatment of non-QCD state 

Lattice QCD provides a powerful tool to calculate the matrix elements with hadronic ini- 
tial/final state. By studying the Euclidean time dependence of the correlation function, we are 
able to pick up the hadronic state of interest. However, this method does not work for the matrix 
elements with non-QCD state. Take the photon state as an example. By using the interpolating 
operator with quantum number J PC = 1 , we expect to extract vector-meson state rather than 
the photon state from the correlation function. To address this problem, Ji and Jung proposed 
an analytic continuation method, which treats the photon as a superposition of a complete set of 
hadron states with appropriate quantum numbers [10]. In our calculation the key observable is the 
matrix element {Y(pii^i)Y(P2i^a)\^ (#)) with a two-photon final state, for which we apply this 
technique. 

We follow the procedure given in Refs. [11, 12]. Using the LSZ reduction formula, we express 
the matrix element in terms of the time-ordered correlation function 

(r(Phh)r(P2,h)\K°(q)) = - lira e*Jp u Xi)e*{p 2 ,k 2 ) 

xpfp'ij d A xd A ye i ^ +i ^{0\T{A^{x)A v {y)}\n\q)) . (3.1) 

At the leading order of perturbative QED expansion, the photon field in the interaction term Hi nt = 
e J d 4 xA^{x)j tl {x) can be contracted with these photon fields existing in Eq. (3.1). After the Wick 
contraction we have 

(r(Puh)r(p2,h)\n°(q)}= -e 2 lim eUpuXiKipiM) 

Pip-tpia 

x pfpf J d 4 x d 4 y d\ d 4 w e^+'i^D^ (x, z)D va (y, w){Q.\T{ j p (z) j a (w) } 1 71° (q) ) .(3.2) 

In Eq. (3.2) j p and j a are the hadronic components of the electromagnetic vector current. The 
photon propagator D^ p (x,z) = —ig >ip J (|^4 e ^ ie ' cancels the inverse propagators outside the 
integral. We then have 

{Y(PiAi)Y(P2,h)\n (q)) = -ie 2 e*(p h X l )e;(p2,X 2 )M tlv ip 1 ,p 2 ) , 
M^ v (p uP2 ) = ij d 4 xe^ x (n\T{j^x)j v (0)}\n Q (q)) . (3.3) 
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where M tlv {p\,p2) is a hadronic matrix element. The form factor ^^yy(rn^,p\,p2) is defined in 
terms of M ilv (p\ ,p2) as 

^ 7l o ry (ml,p 2 l ,pl)=M^ v (p u p 2 )/£^ V apPiP2 > (3-4) 

where the factor £^ va pPiP2 * s induced by the negative parity of K°. 

By an analytic continuation of (3.3) from the Minkowski to Euclidean space-time, we write 

M^( Pl ,p 2 )= lim l — fdhe^-^C^ih,^), (3.5) 

C^h t h t tn) = f d 3 xe-^ f J 3 z^ z (a|r{^(x, ?1 ) 7v (0,f 2 )7r (?,^)}|n) , (3.6) 

where C^ v (?i,^2j^) is a correlation function defined in Euclidean space-time and thus can be cal- 
culated using lattice QCD. The operator J d 3 z e' q ' z K {z,t n ) produces a neutral pion with a spatial 
momentum q. Its amplitude and energy in the ground state are denoted by (j> n ^ and E n ^, respec- 
tively. The four-momentum of the first photon p\ = (co,p\) is chosen as input, while the momen- 
tum of the second photon is given as p% = (E x a— (0,q — p\) by momentum conservation. When 
the squared momentum p \ or p\ exceeds the hadron production threshold, the photon state mixes 
with these hadron states and the analytic continuation fails. To avoid this, we restrict the photon 
momentum in the region p\ 2 < My, where My is the invariant mass of the lowest energy state in 
the vector channel. Although e®^* 2 ' becomes infinitely large when t\ — 1 2 — > °°, a suppression 
factor e -£ v(fi-fc) f r0 m C^, v (ti,t2,t K ) makes the integratal (3.5) convergent. 



4. Lattice setup 

In this calculation we use 2+1 -flavor overlap fermion configurations generated by the JLQCD 
and TWQCD Collaborations [13, 14]. Using the overlap fermions ensures the exact chiral sym- 
metry at even finite lattice spacings. We use a sequence of ensembles with a lattice spacing of 
a = 0. 1 1 fm. The pion mass ranges from 290 to 540 MeV with degenerate up and down quarks. 
The strange quark mass is fixed to be very close to its estimated physical value. The lattice size 
is L 3 x T /a 4 = 16 3 x 48. At two smallest pion masses we also use a larger lattice size L/a = 24 
to check the finite-size (FS) effects. The gauge ensembles are generated by fixing the (global) 
topological charge Q, which results in a finite volume effect of 0(l/L 3 T) [15]. We check the 
significance of this effect by comparing the results with two different values Q = 0, 1. 

We use the all-to-all propagator to construct the correlation function. Since the propagator con- 
tains the information from any source site to any sink site, we are allowed to calculate C^ v {h , hJn) 
at any time slices of t\, ti and t n . Besides, we are able to compute the disconnected diagrams with- 
out extra computer resources. The electromagnetic current is implemented on the lattice as a 
local operator with a renormalization factor calculated nonperturbatively in [16]. 



5. Analysis 

From the large t\ 2 ~ in behavior of C^vihihitx), it is possible to extract the 7T°-ground state. 
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Figure 1: The amplitude A n (z) as a function of T for Figure 2: Left: Contour of (.pf,/>§) rescaled 
momentum setups 1 (left) and 2 (right). The black by 1 /My for our momentum setups. Right: 
(red) curves indicate the lattice (VMD) amplitudes. F^o^m^^p^p^) as a function of p\jM v . 



We define the amplitude A n as 

A»(t)= lim C^h^t^/e-^-^ , (5.1) 
with T = ?i — ?2 and ? = minjfi^}, and obtain M /lv (pi,P2) by performing an integral 

^ f jfdT + 1° ^-^^(t)) . (5.2) 

The spatial momenta are assigned as p\ = ^(0,0,0), q = ^(0,0, 1) (setup 1) and p\ = ^(0,0, 1), 
q = ^ (0,0,0) (setup 2). The resulting amplitudes A 7[ (z) for these setups are shown by the black 
curves in Fig. 1. The amplitude A^ md (t), constructed from the vector-meson-dominance (VMD) 
form factor i^K,/, 2 ,/? 2 ) = c v Gv{p\)Gv(p\), with G v (p 2 ) = M 2 /{M 2 - p 2 ) the vector me- 
son propagator and cy a constant is also plotted by red curves. We find that A^ md (t) give a good 
approximation to the lattice data at large |t| but fails to match them at small |t|. This is because no 
information of the vector- meson excited states are contained in ^^o^( m |>/ 7 i>/ 7 2)- Although the 
lattice data are truncated due to the finite time extent T, we are still able to perform the integral (5.2) 
from — oo to +oo since the A k {t) at the large |t| is dominated by the lowest vector meson. 

When performing the integral (5.2) we can tune the value of the photon energy ft) continu- 
ously. As shown in the left panel of Fig. 2, a pair {p\,p\) = (ft) 2 — p 2 , (E n ^ — ft)) 2 — (q — pi) 2 ) 
forms a continuous contour on the {p\,p\) plane for p\ 2 < My/2. Evaluating ^^^(m^, p 2 , p 2 ) 
along this contour, we obtain the data plotted in the right panel of Fig. 2. We use the fit ansatz 
^ n {ml,p\,pl) = 

c v G v (p 2 )G v (p 2 )+l,c m {(p 2 ) m G v (p 2 ) + (p 2 ) m G v (pty (5.3) 

m m.n 

which includes the contribution from the lowest vector meson through Gv{p 2 2 ) and the contri- 
bution from excited states as a polynomial of p\ 2 . We perform the combined fit of the lattice 
data to Eq. (5.3) with four free parameters: cy, Co, co.o and co,i = ci i0 . The contributions from 
the higher-order terms are not significant when compared to the statistical error and thus can be 
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Figure 3: F(n%,0,Q), 

8V> gvny, F n and FS corrected F(m|,0,0) as a function of nv\ from top to bottom 
panels. In each panel, data with (L/a,Q) = (16,0), (24,0) and (16,1) are plotted by the blue, red and 
green symbols, respectively. The yellow symbols indicate the Particle Data Group (PDG) [17] or PrimEx [6] 
experimental values for reference. The solid (dashed) curves show the result of the fit to the linear (quadratic) 
function. The dataset used in the fit is explained in the text. 

neglected. The fitting curves are shown in the right panel of Fig. 2. We find that the single for- 
mula (5.3) describes the data with different momentum setups. Using the resulting fit parameters 
and extrapolating ^ n Qyy{m 2 n ,p\,p\) to the soft photon limit, we obtain the normalized form factors 
F(m\, 0,0) = {A% 1 F n )^ 7 fiyy(jn 1 n ^^\ which are plotted in the uppermost panel of Fig. 3. 

Using the results of F(m|,0,0) we analyze the various systematic effects including FS effects, 
fixed-topology effects, disconnected-diagram effects, higher-order effects in chiral extrapolation 
and possible discretization effects in the evaluation of Eq. (3.5). For more details of the analysis, 
we refer the readers to our recent publication [18]. 

Among the various systematic effects, the largest one is the conventional FS effect. As shown 
in Fig. 3, at m% « 290 MeV we find that F(m 2 , 0,0) calculated at L/a = 16 lattice is 27% less 
than the one at L/a = 24. To qualitatively understand this FS effect, we insert the ground state 
into (j^ijv^ ) and approximate this three-point correlation function with three hadronic matrix 
elements: (j^,j v 7l ; ) — > (Q.\j^\V,e)(V,e\j v \7i )(^°|^°|Q). These matrix elements are related to 
the electromagnetic coupling gy, the Vny coupling gy Ky and the pion decay constant F K . In our 
calculation we do not observe significant FS effect in My but find 8%, 7% and 9% shifts in gy, gy K y 
and F K , respectively, from L/a = 16 to 24, as shown in Fig. 3. These FS effects may accumulate in 
the three-point function. We estimate the FS corrections to gy, gy n y and F % by adding a correction 
term, e~ m " L , into the linear fit form in the chiral extrapolation of each quantity or using NNLO 
chiral perturbation theory. We then combine the FS corrections to each quantity as an estimation 
of the total FS correction to F(ra|,0,0). As shown in the lowest panel of Fig. 3 the FS corrected 
F(m|,0,0) alL/a = 16 agrees with those dlL/a = 24. 

We use two methods of treating the FS effects: 1. We use the uncorrected F(m|, 0,0) with 
m n L > 4 to perform the chiral extrapolation. Namely, we exclude the L/a = 16 data points at two 
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lowest pion masses. A linear function in m\ is used as an fit ansatz. 2. We use the FS corrected 
data of all ensembles to perform a linear extrapolation. The difference between the results from the 
two methods is considered as a systematic error. 

We quote our results for F (0,0,0) and T n 0yy in the isospin symmetric limit as F (0,0,0) = 
1.009(22)(29) and T^Oyy = 7.83(31)(49) eV, where the first error is statistical and the second one 
is systematic. Our results reproduce the predication of the ABJ anomaly F (0,0,0) = 1 and agree 
with the PrimEx measurement T^oyy = 7.82(22) eV [6]. For future improvements, isospin breaking 
effects due to the light quark mass difference need to be included. 

Numerical simulations are performed on the Hitachi SRI 6000 at Yukawa Institute of Theoret- 
ical Physics and at High Energy Accelerator Research Organization under a support of its Large 
Scale Simulation Program (No. 11-05). This work is supported in part by the Grant-in-Aid of 
the Japanese Ministry of Education (No.21674002, 21684013), the Grant-in-Aid for Scientific Re- 
search on Innovative Areas (No. 2004: 20105001, 20105002, 20105003, 20105005, 23105710), 
and SPIRE (Strategic Program for Innovative Research). 
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